home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / minix / libsrc~1.z / libsrc~1 / atof.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-12-28  |  12.4 KB  |  543 lines
  1. /*
  2.  *     double strtod (str, endptr);
  3.  *     _CONST char *str;
  4.  *     _CONST char **endptr;
  5.  *        if !NULL, on return, points to char in str where conv. stopped
  6.  *
  7.  *     double atof (str)
  8.  *     _CONST char *str;
  9.  *
  10.  * recognizes:
  11.  *   [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int]]]]
  12.  *
  13.  * returns:
  14.  *    the number
  15.  *        on overflow: HUGE and errno = ERANGE
  16.  *        on underflow: -HUGE and errno = ERANGE
  17.  *
  18.  * bugs:
  19.  *   naive over/underflow detection
  20.  *
  21.  *    ++jrb bammi@dsrgsun.ces.cwru.edu
  22.  *
  23.  * 07/29/89:
  24.  *    ok, before you beat me over the head, here is a hopefully
  25.  *    much improved one, with proper regard for rounding/precision etc
  26.  *    (had to read up on everything i did'nt want to know in K. V2)
  27.  *    the old naive coding is at the end bracketed by ifdef __OLD__.
  28.  *    modeled after peter housels posting on comp.os.minix.
  29.  *    thanks peter!
  30.  *        ++jrb
  31.  */
  32. #if !(defined(unix) || defined(minix))
  33. #include <stddef.h>
  34. #include <stdlib.h>
  35. #include <float.h>
  36. #endif
  37. #include <errno.h>
  38. #include <assert.h>
  39.  
  40. #ifdef minix
  41. #include "lib.h"
  42. #endif
  43.  
  44. extern int errno;
  45. #if (defined(unix))
  46. #define NULL     ((void *)0)
  47. #endif
  48.  
  49. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  50. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  51. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  52. #define Issign(c)    ((c == '-') || (c == '+'))
  53. #define Val(c)        ((c - '0'))
  54.  
  55. #define MAXDOUBLE    DBL_MAX
  56. #define MINDOUBLE    DBL_MIN
  57.  
  58. #define MAXF  1.797693134862316
  59. #define MINF  2.225073858507201
  60. #define MAXE  308
  61. #define MINE  (-308)
  62.  
  63. /* some doubt about this */
  64. static unsigned long __huge[2] = { 0x7ff00000L, 0L }; /* ieee infinity */
  65. #define HUGE  (*((double *)&__huge[0]))
  66.  
  67. /* another alias for ieee double */
  68. struct ldouble {
  69.     unsigned long hi, lo;
  70. };
  71.  
  72. /*
  73.  * mul 64 bit accumulator by 10 and add digit
  74.  */
  75. static int __ten_mul(acc, digit)
  76. double *acc;
  77. int digit;
  78. {
  79.     register unsigned long d0, d1, d2, d3;
  80.     register          short i;
  81.     
  82.     d2 = d0 = ((struct ldouble *)acc)->hi;
  83.     d3 = d1 = ((struct ldouble *)acc)->lo;
  84.  
  85.     /* check possibility of overflow */
  86. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  87.     if( (d0 & 0x70000000L) != 0 )
  88.     /* report overflow somehow */
  89.     return 1;
  90.     
  91.     /* 10acc == 2(4acc + acc) */
  92.     for(i = 0; i < 2; i++)
  93.     {  /* 4acc == ((acc) << 2) */
  94.     asm volatile("    lsll    #1,%1;
  95.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  96.         : "=d" (d0), "=d" (d1)
  97.         : "0"  (d0), "1"  (d1) );
  98.     }
  99.  
  100.     /* 4acc + acc */
  101.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  102.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  103.  
  104.     /* (4acc + acc) << 1 */
  105.     asm volatile("  lsll    #1,%1;
  106.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  107.         : "=d" (d0), "=d" (d1)
  108.         : "0"  (d0), "1"  (d1) );
  109.  
  110.     /* add in digit */
  111.     d2 = 0;
  112.     d3 = digit;
  113.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  114.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  115.  
  116.  
  117.     /* stuff result back into acc */
  118.     ((struct ldouble *)acc)->hi = d0;
  119.     ((struct ldouble *)acc)->lo = d1;
  120.  
  121.     return 0;    /* no overflow */
  122. }
  123.  
  124. #include "flonum.h"
  125.  
  126. double __adjust(acc, dexp, sign)
  127. double *acc;    /* the 64 bit accumulator */
  128. int     dexp;    /* decimal exponent       */
  129. int    sign;    /* sign flag          */
  130. {
  131.     register unsigned long  d0, d1, d2, d3;
  132.     register          short i;
  133.     register           short bexp = 0; /* binary exponent */
  134.     unsigned short tmp[4];
  135.     double r;
  136.     struct bitdouble *rp = (struct bitdouble *)&r;
  137. #ifdef __STDC__
  138.     double norm( double d, int exp, int sign, int rbits );
  139. #else
  140.     extern double norm();
  141. #endif    
  142.     d0 = ((struct ldouble *)acc)->hi;
  143.     d1 = ((struct ldouble *)acc)->lo;
  144.     while(dexp != 0)
  145.     {    /* something to do */
  146.     if(dexp > 0)
  147.     { /* need to scale up by mul */
  148.         while(d0 > 429496729 ) /* 2**31 / 5 */
  149.         {    /* possibility of overflow, div/2 */
  150.         asm volatile(" lsrl    #1,%1;
  151.                     roxrl    #1,%0"
  152.             : "=d" (d1), "=d" (d0)
  153.             : "0"  (d1), "1"  (d0));
  154.         bexp++;
  155.         }
  156.         /* acc * 10 == 2(4acc + acc) */
  157.         d2 = d0;
  158.         d3 = d1;
  159.         for(i = 0; i < 2; i++)
  160.         {  /* 4acc == ((acc) << 2) */
  161.         asm volatile("    lsll    #1,%1;
  162.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  163.                  : "=d" (d0), "=d" (d1)
  164.                  : "0"  (d0), "1"  (d1) );
  165.         }
  166.  
  167.         /* 4acc + acc */
  168.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  169.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  170.  
  171.         /* (4acc + acc) << 1 */
  172.         bexp++;    /* increment exponent to effectively acc * 10 */
  173.         dexp--;
  174.     }
  175.     else /* (dexp < 0) */
  176.     { /* scale down by 10 */
  177.         while((d0 & 0xC0000000L) == 0)
  178.         {    /* preserve prec by ensuring upper bits are set before div */
  179.         asm volatile(" lsll  #1,%1;
  180.                     roxll #1,%0" /* shift L to move bits up */
  181.             : "=d" (d0), "=d" (d1)
  182.             : "0"  (d0), "1"  (d1) );
  183.         bexp--;    /* compensate for the shift */
  184.         }
  185.         /* acc/10 = acc/5/2 */
  186.         *((unsigned long *)&tmp[0]) = d0;
  187.         *((unsigned long *)&tmp[2]) = d1;
  188.         d2 = (unsigned long)tmp[0];
  189.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  190.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  191.         for(i = 1; i < 4; i++)
  192.         {
  193.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  194.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  195.         tmp[i] = (unsigned short)d2;
  196.         }
  197.         d0 = *((unsigned long *)&tmp[0]);
  198.         d1 = *((unsigned long *)&tmp[2]);
  199.         /* acc/2 */
  200.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  201.         dexp++;
  202.     }
  203.     }
  204.     
  205.     /* stuff the result back into acc */
  206.     ((struct ldouble *)acc)->hi = d0;
  207.     ((struct ldouble *)acc)->lo = d1;
  208.  
  209.     /* normalize it */
  210.     r = norm( *acc, ((0x03ff - 1) + 53), sign, 0 );
  211.     /* now shove in the binary exponent */
  212.     rp->exp += bexp;
  213.     return r;
  214. }
  215.  
  216. /* flags */
  217. #define SIGN    0x01
  218. #define ESIGN    0x02
  219. #define DECP    0x04
  220.  
  221. double strtod (s, endptr)
  222. _CONST register char *s;
  223. _CONST char **endptr;
  224. {
  225.     double         accum = 0.0;
  226.     register short flags = 0;
  227.     register short texp  = 0;
  228.     register short e     = 0;
  229. #ifdef __STDC__
  230.     int __ten_mul(double *, int);
  231.     double __adjust(double *, int, int);
  232. #else
  233.     extern int __ten_mul();
  234.     extern double __adjust();
  235. #endif
  236.     
  237.     assert ((s != NULL));
  238.     errno = 0;
  239.     
  240.     while(Isspace(*s)) s++;
  241.     if(*s == '\0')
  242.     {    /* just leading spaces */
  243.     if(endptr != NULL) *endptr = s;
  244.     return 0.0;
  245.     }
  246.     
  247.     
  248.     if(Issign(*s))
  249.     {
  250.     if(*s == '-') flags = SIGN;
  251.     if(*++s == '\0')
  252.     {   /* "+|-" : should be an error ? */
  253.         if(endptr != NULL) *endptr = s;
  254.         return 0.0;
  255.     }
  256.     }
  257.     
  258.     for(; (Isdigit(*s) || (*s == '.')); s++)
  259.     {
  260.     if(*s == '.')
  261.         flags |= DECP;
  262.     else
  263.     {
  264.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  265.         if(flags & DECP) texp--;
  266.     }
  267.     }
  268.  
  269.     if(Ise(*s))
  270.     {
  271.     if(*++s != '\0') /* skip e|E|d|D */
  272.     {  /* ! ([s]xxx[.[yyy]]e)  */
  273.     
  274.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  275.         if(*s != '\0')
  276.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  277.     
  278.         if(Issign(*s))
  279.             if(*s++ == '-') flags |= ESIGN;
  280.     
  281.         if(*s != '\0')
  282.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  283.  
  284.             for(; Isdigit(*s); s++)
  285.             e = (((e << 2) + e) << 1) + Val(*s);
  286.  
  287.             /* dont care what comes after this */
  288.             if(flags & ESIGN)
  289.             texp -= e;
  290.             else
  291.             texp += e;
  292.         }
  293.         }
  294.     }
  295.     }
  296.  
  297.     if(endptr != NULL) *endptr = s;
  298.     if(accum == 0.0) return 0.0;
  299.     
  300.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  301.     
  302. }
  303.  
  304. double atof(s)
  305. _CONST char *s;
  306. {
  307.     return strtod(s, (_CONST char **)NULL);
  308. }
  309.  
  310. #ifdef TEST
  311. #ifdef __MSHORT__
  312. #error "please run this test in 32 bit int mode"
  313. #endif
  314.  
  315. #define NTEST 10000L
  316.  
  317. #ifdef unix
  318. #ifdef __MSHORT__
  319. #define    RAND_MAX    (0x7FFF)    /* maximum value from rand() */
  320. #else
  321. #define    RAND_MAX    (0x7FFFFFFFL)    /* maximum value from rand() */
  322. #endif
  323. #endif
  324.  
  325. main()
  326. {
  327.     
  328.     double expected, result, e, max_abs_err;
  329.     char buf[128];
  330.     register long i, errs;
  331.     register int s;
  332. #ifdef __STDC__
  333.     double atof(_CONST char *);
  334.     int rand(void);
  335. #else
  336.     extern double atof();
  337.     extern int rand();
  338. #endif
  339.  
  340. #if 0
  341.     expected = atof("3.14159265358979e23");
  342.     expected = atof("3.141");
  343.     expected = atof(".31415"); 
  344.     printf("%f\n\n", expected);
  345.     expected = atof("3.1415"); 
  346.     printf("%f\n\n", expected);
  347.     expected = atof("31.415"); 
  348.     printf("%f\n\n", expected);
  349.     expected = atof("314.15"); 
  350.     printf("%f\n\n", expected);
  351.  
  352.     expected = atof(".31415"); 
  353.     printf("%f\n\n", expected);
  354.     expected = atof(".031415"); 
  355.     printf("%f\n\n", expected);
  356.     expected = atof(".0031415"); 
  357.     printf("%f\n\n", expected);
  358.     expected = atof(".00031415"); 
  359.     printf("%f\n\n", expected);
  360.     expected = atof(".000031415"); 
  361.     printf("%f\n\n", expected);
  362.  
  363.     expected = atof("-3.1415e-9"); 
  364.     printf("%20.15e\n\n", expected);
  365.  
  366.     expected = atof("+3.1415e+009"); 
  367.     printf("%20.15e\n\n", expected);
  368. #endif
  369.  
  370.     expected = atof("+3.123456789123456789"); 
  371.     printf("%30.25e\n\n", expected);
  372.  
  373.     expected = atof(".000003123456789123456789"); 
  374.     printf("%30.25e\n\n", expected);
  375.  
  376.     expected = atof("3.1234567891234567890000000000"); 
  377.     printf("%30.25e\n\n", expected);
  378.  
  379. #if 0
  380.     expected = atof("1.7e+308");
  381.     if(errno != 0)
  382.     {
  383.     printf("%d\n", errno);
  384.     }
  385.     else    printf("1.7e308 OK %g\n", expected);
  386.     expected = atof("1.797693e308");    /* anything gt looses */
  387.     if(errno != 0)
  388.     {
  389.     printf("%d\n", errno);
  390.     }
  391.     else    printf("Max OK %g\n", expected);
  392.     expected = atof("2.225073858507201E-307");
  393.     if(errno != 0)
  394.     {
  395.     printf("%d\n", errno, expected);
  396.     }
  397.     else    printf("Min OK %g\n", expected);
  398. #endif
  399.     
  400.     max_abs_err = 0.0;
  401.     for(errs = 0, i = 0; i < NTEST; i++)
  402.     {
  403.     expected = (double)(s = rand()) / (double)rand();
  404.     if(s > (RAND_MAX >> 1)) expected = -expected;
  405.     sprintf(buf, "%.14e", expected);
  406.     result = atof(buf);
  407.     e = (expected == 0.0) ? result : (result - expected)/expected;
  408.     if(e < 0) e = (-e);
  409.     if(e > 1.0e-6) 
  410.     {
  411.         errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
  412.     }
  413.     if (e > max_abs_err) max_abs_err = e;
  414.     }
  415.     printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
  416. }
  417. #endif /* TEST */
  418.  
  419. /* old naive coding */
  420. #ifdef __OLD__
  421. static double __ten_pow(r, e)
  422. double r;
  423. register int e;
  424. {
  425.     if(e < 0)
  426.     for(; e < 0; e++) r /= 10.0;
  427.     else
  428.     for(; e > 0; --e) r *= 10.0;
  429.     return r;
  430. }
  431.  
  432. #define RET(X)     {goto ret;}
  433.  
  434. double strtod (s, endptr)
  435. _CONST register char *s;
  436. _CONST char **endptr;
  437. {
  438.     double f = 0.1, r = 0.0, accum = 0.0;
  439.     register int  e = 0, esign = 0, sign = 0;
  440.     register int texp = 0, countexp;
  441. #ifdef __STDC__
  442.     double __ten_pow(double, int);
  443. #else
  444.     extern double __ten_pow();
  445. #endif
  446.     
  447.     assert ((s != NULL));
  448.     errno = 0;
  449.     
  450.     while(Isspace(*s)) s++;
  451.     if(*s == '\0') RET(r);    /* just spaces */
  452.     
  453.     if(Issign(*s))
  454.     {
  455.     if(*s == '-') sign = 1;
  456.     s++;
  457.     if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
  458.     }
  459.     countexp = 0;
  460.     while(Isdigit(*s))
  461.     {
  462.     if(!countexp && (*s != '0')) countexp = 1;
  463.     accum = (accum * 10.0) + Val(*s);
  464.     /* should check for overflow here somehow */
  465.     s++; 
  466.     if(countexp) texp++;
  467.     }
  468.     r = (sign ? (-accum) : accum);
  469.     if(*s == '\0') RET(r);  /* [s]xxx */
  470.     
  471.     countexp = (texp == 0);
  472.     if(*s == '.')
  473.     {
  474.     s++;
  475.     if(*s == '\0') RET(r); /* [s]xxx. */
  476.     if(!Ise(*s))
  477.     {
  478.         while(Isdigit(*s))
  479.         {
  480.         if(countexp && (*s == '0')) --texp;
  481.         if(countexp && (*s != '0')) countexp = 0;
  482.         accum = accum + (Val(*s) * f);
  483.         f = f / 10.0;
  484.         /* overflow (w + f) ? */
  485.         s++;
  486.         }
  487.         r = (sign ? (-accum) : accum);
  488.         if(*s == '\0') RET(r); /* [s]xxx.yyy  */
  489.     }
  490.     }
  491.     if(!Ise(*s)) RET(r);    /* [s]xxx[.[yyy]]  */
  492.     
  493.     s++; /* skip e */
  494.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e  */
  495.     
  496.     while(Isspace(*s)) s++;
  497.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space]  */
  498.     
  499.     if(Issign(*s))
  500.     {
  501.     if(*s == '-') esign = 1;
  502.     s++;
  503.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s]  */
  504.     }
  505.     
  506.     while(Isdigit(*s))
  507.     {
  508.     e = (e * 10) + Val(*s);
  509.     s++;
  510.     }
  511.     /* dont care what comes after this */
  512.     if(esign) e = (-e);
  513.     texp += e;
  514.     
  515.     /* overflow ? */ /* FIXME */
  516.     if( texp > MAXE)
  517.     {
  518.     if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
  519.     {
  520.         errno = ERANGE;
  521.         r = ((sign) ? -HUGE : HUGE);
  522.         RET(r);
  523.     }
  524.     }
  525.     
  526.     /* underflow -- in reality occurs before this */ /* FIXME */
  527.     if(texp < MINE)
  528.     {
  529.     errno = ERANGE;
  530.     r = ((sign) ? -HUGE : HUGE);
  531.     RET(r);
  532.     }
  533.     r = __ten_pow(r, e);
  534.     /* fall through */
  535.     
  536.     /* all returns end up here, with return value in r */
  537.   ret:
  538.     if(endptr != NULL)
  539.     *endptr = s;
  540.     return r;
  541. }
  542. #endif /* __OLD__ */
  543.